↜ Back to index Introduction to Numerical Analysis 2

In the previous lecture, we started to speed up the conjugate gradient method for sparse matrices.

A matrix is said to be sparse (疎行列) if most of its entries are zero. If a matrix A is sparse, we have an opportunity to save memory necessary to store A since we need to only remember the nonzero entries. Furthermore, we might be able to speed up the computation of the matrix-vector product Ax, (Ax)_i = \sum_{j=1}^n A_{ij}x_j, by not performing the multiplication and addition on the right-hand side whenever A_{ij} is zero.

As we saw in the implementation of lap_mul in Lecture 5, if the distribution of the nonzero elements in a matrix A is regular, we do not even need to remember the matrix in a variable! But of course life is not always that easy.

Let us consider a sparse matrix A = \begin{pmatrix} 7 & 0 & -3 & 0\\ 2 & 0 & 0 & 0\\ 0 & 4 & 0 & 0\\ 0 & 0 & -2 & 0\\ \end{pmatrix} This matrix has 5 nonzero entries out of the total of 16 entries. Let us consider how we can take an advantage of this fact to speed up the computation of Ax and save some memory.

Our goal is to store the nonzero entries in such a way that we can easily compute the matrix-vector product Ax given as (Ax)_i = \sum_{j=1}^n A_{ij} x_j. Whenever A_{ij} is zero for some i and j we can skip that term. So we need to remember only the values of nonzero entries and the row and column in which they are located. This is enough to compute Ax.

The most straightforward way is to use 3 arrays of size 5:

integer, parameter :: k = 5
real v(k) 
integer ri(k), ci(k)

We use v to store the values of the nonzero elements, and ri and ci to store the row and column indices. For example, we can set the values to

 v = [7., -3., 2., 4., -2.]
 ri = [1, 1, 2, 3, 4]
 ci = [1, 3, 1, 2, 3]

This way we only need to remember 3 numbers for each nonzero element of A. This type of sparse matrix storage is called coordinate list.

Let us now try to implement the matrix-vector multiplication Ax given the 3 arrays above.

One implementation, assuming that A is a square matrix, can look like:

function sparse_coo_mul(v, ri, ci, x) result(y)
    real v(:), x(:), y(size(x))
    integer ri(:), ci(:)
    integer i

    y = 0
    do i = 1,size(v)
        y(ri(i)) = y(ri(i)) + v(i) * x(ci(i))
    enddo
end

Note that we need to pass all 3 arrays v, ri and ci to specify the matrix.

We can further improve the storage of a sparse matrix. Here we consider the nonzero elements to be ordered according to their row index. This way the content of the array ri for matrix M above is

ri = [1, 1, 3, 3, 4, 4, 4, 5]

Note that many of the numbers are repeated. To save space, we can instead just remember which entries belong to the first row, to the second row, and so on. An easy way is to replace ri by an array rb of size n + 1, where n is the number of rows, and set rb(i) to be the index of the first nonzero entry belonging to row i. If a row i is empty, we set rb(i) to be equal to rb(i+1). We set rb(n+1) to k + 1. For matrix M we get

rb = [1, 3, 3, 5, 8, 9]

The number of nonzero entries in row i is given as rb(i+1) - rb(i). Note that both rb(2) and rb(3) are 3. This indicates that row 2 of the matrix has no nonzero entries.

To summarize, this way we can represent the matrix M using the following 3 arrays.

 integer, parameter :: k = 8
 integer, parameter :: n = 5
 real v(k) 
 integer rb(n + 1), ci(k)

 v = [5., 9., -2., -7., 6., 3., -8., 2.]
 rb = [1, 3, 3, 5, 8, 9]
 ci = [2, 3, 1, 4, 2, 3, 4, 1]

This type of sparse matrix storage is called compressed sparse row.

Once again, we are interested in the matrix-vector product. Using this format is simple.

function sparse_csr_mul(v, rb, ci, x) result(y)
    real v(:), x(:), y(size(x))
    integer rb(:), ci(:)
    integer i,j

    do i = 1,size(x)
        y(i) = 0
        do j = rb(i),rb(i+1) - 1
            y(i) = y(i) + v(j) * x(ci(j)) 
        enddo
    enddo
end

Check for yourself that the above code computes the correct product for the example matrix above.

Let us now consider examples of how to initialize the sparse matrix storage for a practical matrix whose size depends on the input value.

implicit none
real, dimension(:), allocatable :: v
integer, dimension(:), allocatable :: rb, ci
integer i, j, n, k, t

read *, n

! number of nonzero elements
k = 3 * n - 2

allocate(v(k))
allocate(rb(n + 1))
allocate(ci(k))

! t is the index where the next element will be inserted
t = 1
do i = 1, n
    rb(i) = t

    if (i > 1) then
        ci(t) = i - 1
        v(t) = -1
        t = t + 1
    endif

    ci(t) = i
    v(t) = 2.
    t = t + 1

    if (i < n) then
        ci(t) = i + 1
        v(t) = -1
        t = t + 1
    endif
enddo
rb(n + 1) = k + 1

print *, 'v', v
print *, 'rb', rb
print *, 'ci', ci
end

$ ./a.out <<< "1"
 v   2.00000000    
 rb           1           2
 ci           1
$ ./a.out <<< "2"
 v   2.00000000      -1.00000000      -1.00000000       2.00000000    
 rb           1           3           5
 ci           1           2           1           2
$ ./a.out <<< "3"
 v   2.00000000      -1.00000000      -1.00000000       2.00000000      -1.00000000      -1.00000000       2.00000000    
 rb           1           3           6           8
 ci           1           2           1           2           3           2           3

Exercises

$ ./a.out <<< "2"
 v   1.00000000       1.00000000       1.00000000       1.00000000    
 rb           1           3           5
 ci           1           2           1           2

$ ./a.out <<< "3"
 v   1.00000000       1.00000000       1.00000000       1.00000000       1.00000000    
 rb           1           3           4           6
 ci           1           3           2           1           3

$ ./a.out <<< "4"
 v   1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000       1.00000000    
 rb           1           3           5           7           9
 ci           1           4           2           3           2           3           1           4